## Setup ##

rm(list=ls())
setwd("---")#set working directory, and create a sub-folder called "Figures"

# Packages

ipak <- function(pkg){new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}

packages <- c("tidyverse", "tidytext", "lubridate", "stringr",
              "kableExtra", "scales", "janitor","ggplot2","ggstance","descr")

ipak(packages)


# Load in and clean datasets

civ <- read.csv("civ_survey.csv",stringsAsFactors = F)
civ$complete <- !(is.na(civ$vig2_dem_army)&is.na(civ$vig2_gop_army))
civ <- civ[civ$complete==T,] #remove incomplete observations

civ <- civ %>% 
  janitor::clean_names() %>% 
  mutate_at(vars(age:vig2_gop_militia), ~na_if(., 99))  %>% # 99 did not answer
  mutate_at(vars(age:vig2_gop_militia), ~na_if(., 667)) # typo in age category

mil <- read.csv("mil_survey.csv",stringsAsFactors = F)
mil$complete <- !(is.na(mil$vig1_dem_army)&is.na(mil$vig1_gop_army))
mil <- mil[mil$complete==T,] # remove incomplete observations

# merge datasets

mil$military <- 1
civ$military <- 0
comp_dat <- merge(civ,mil,all = T,by=intersect(names(civ),names(mil)))


# Variables

# PID scale (1=strong dem, 7=strong rep)
comp_dat$pid_scale <- NA
comp_dat$pid_scale[!is.na(comp_dat$pid)&comp_dat$pid==2] <- comp_dat$dem_pt_2[!is.na(comp_dat$pid)&comp_dat$pid==2]
comp_dat$pid_scale[!is.na(comp_dat$pid)&comp_dat$pid==3] <- ifelse(comp_dat$indepen_pt_2[!is.na(comp_dat$pid)&comp_dat$pid==3]==1,5,
                                              ifelse(comp_dat$indepen_pt_2[!is.na(comp_dat$pid)&comp_dat$pid==3]==2,3,4))
comp_dat$pid_scale[!is.na(comp_dat$pid)&comp_dat$pid==1] <- 8-comp_dat$gop_pt_2[!is.na(comp_dat$pid)&comp_dat$pid==1]

# PID binary w/ Leaners (0=Dems, 1=Reps, NA=other)
comp_dat$pid_lean <- ifelse(!is.na(comp_dat$pid)&comp_dat$pid_scale<=3,0,
                           ifelse(!is.na(comp_dat$pid)&comp_dat$pid_scale>=5,1,NA))

# PID binary w/o Leaners (1=Dems, 2=Indep, 3=Reps, NA=other)
comp_dat$pid_bin <- ifelse(!is.na(comp_dat$pid)&comp_dat$pid==2,1,
                            ifelse(!is.na(comp_dat$pid)&comp_dat$pid==3,2,
                                   ifelse(!is.na(comp_dat$pid)&comp_dat$pid==1,3,NA)))

# vignette responses

# treatments (dem protesters = 0, gop protesters = 1)
comp_dat$treat_1 <- ifelse(!is.na(comp_dat$vig1_dem_army),1,0)
comp_dat$treat_2 <- ifelse(!is.na(comp_dat$vig2_dem_army),0,1)

# treatments (same party protesters=0,other=1)
comp_dat$same_1 <- ifelse(comp_dat$pid_lean==comp_dat$treat_1,1,0)
comp_dat$same_2 <- ifelse(comp_dat$pid_lean==comp_dat$treat_2,1,0)

#all responses merged
comp_dat$vig1_army <- ifelse(comp_dat$treat_1==1,comp_dat$vig1_dem_army,comp_dat$vig1_gop_army)
comp_dat$vig1_police_eq <- ifelse(comp_dat$treat_1==1,comp_dat$vig1_dem_police_eq,comp_dat$vig1_gop_police_eq)
comp_dat$vig1_militia <- ifelse(comp_dat$treat_1==1,comp_dat$vig1_dem_militia,comp_dat$vig1_gop_militia)
comp_dat$vig2_army <- ifelse(comp_dat$treat_2==0,comp_dat$vig2_dem_army,comp_dat$vig2_gop_army)
comp_dat$vig2_police_eq <- ifelse(comp_dat$treat_2==0,comp_dat$vig2_dem_police_eq,comp_dat$vig2_gop_police_eq)
comp_dat$vig2_militia <- ifelse(comp_dat$treat_2==0,comp_dat$vig2_dem_militia,comp_dat$vig2_gop_militia)


#### Descriptive Stats ####

table(civ$pid)
civ$covariate <- 1
civ$pid[civ$pid==4] <- "None"
civ$gender <- ifelse(is.na(civ$gender),"Non-Binary",ifelse(civ$gender==1, "Female","Male"))
civ$age_grp = case_when(
  civ$age < 35 ~ "< 35",
  civ$age >= 35 & civ$age < 55 ~ "35-55",
  civ$age >= 55 ~ "55 >")
comp_dat$service_branch <- ifelse(comp_dat$military==1,comp_dat$service,
                                  ifelse(is.na(comp_dat$service_branch),"None",comp_dat$service_branch))

# Table A2 (Civilian sample only)

freq(civ$gender)
freq(civ$age_grp)
freq(civ$education)
freq(civ$race)
freq(civ$work)
freq(civ$location)
freq(civ$mil_connections)
freq(civ$vet_status)

# Table A3 (Civilian and Military Sample)

crosstab(comp_dat$pid_bin,comp_dat$military,prop.c = TRUE)
freq(mil$tenure)
crosstab(comp_dat$service_branch,comp_dat$military,prop.c = TRUE)


#### Main Results ####


# Figure A2: graph of treatment effects by party

results <- data.frame(Vignette=numeric(length=12),DV=character(length=12),DV_Name=character(length=12),
                      Respondents=character(length=12),ATE=numeric(length=12),SE=numeric(length=12))

results$Vignette <- rep(1:2, each=6)
results$Respondents <- rep(c("Dems/Lean","GOP/Lean"),times=6)
results$DV_Name <- rep(rep(c("Deploy Army","Militarize Police","Armed Volunteers"),each=2),times=2)
results$DV <- rep(1:6,each=2)

all_vars <- cbind(comp_dat$vig1_army,ifelse(!is.na(comp_dat$vig1_police_eq),ifelse(comp_dat$vig1_police_eq%in%c(1,3),1,0),NA),comp_dat$vig1_militia,
                  comp_dat$vig2_army,ifelse(!is.na(comp_dat$vig1_police_eq),ifelse(comp_dat$vig2_police_eq%in%c(1,3),1,0),NA),comp_dat$vig2_militia)
all_treats <- cbind(comp_dat$treat_1,comp_dat$treat_2)
all_vars[,c(1,3,4,6)] <- (all_vars[,c(1,3,4,6)]-1)/4

covariates <- cbind(comp_dat$age[],comp_dat$gender==1,comp_dat$education,as.factor(comp_dat$race),comp_dat$vet_status,
                    comp_dat$trust_prez,comp_dat$trust_congress,comp_dat$trust_police,comp_dat$trust_military)

for(i in 1:12){
  dv <- all_vars[comp_dat$military==0&comp_dat$pid_lean==ifelse(results$Respondents[i]=="GOP/Lean",1,0),results$DV[i]]
  treat <- all_treats[comp_dat$military==0&comp_dat$pid_lean==ifelse(results$Respondents[i]=="GOP/Lean",1,0),results$Vignette[i]]
  covar <- covariates[comp_dat$military==0&comp_dat$pid_lean==ifelse(results$Respondents[i]=="GOP/Lean",1,0),]
  results[i,5:6] <- summary(lm(dv~treat+covar[,1]+covar[,2]+covar[,3]+covar[,4]+covar[,5]+covar[,6]+covar[,7]+covar[,8]+covar[,9]))$coefficients[2,1:2]
}

results$Upper <- results$ATE+1.96*results$SE
results$Lower <- results$ATE-1.96*results$SE

results$Respondents <- factor(results$Respondents,levels=c("GOP/Lean","Dems/Lean"))
results$DV_Name <- factor(results$DV_Name,levels=c("Armed Volunteers","Militarize Police","Deploy Army"))
results$Vignette_Name <- ifelse(results$Vignette==1,"Capitol Protest","Downtown Riot")

pdf("Figures/Figure_A2.pdf", width = 6, height = 3)
ggplot(results, aes(y = DV_Name, x = ATE, color=Respondents)) +
  geom_vline(xintercept = 0, linetype = 2) + 
  facet_grid(Vignette_Name ~., scale="free_y", space = "free_y") +
  geom_point(size = 2,position=position_dodgev(height=-0.5)) + 
  scale_x_continuous(limits = c(-.2, .2)) +
  geom_errorbarh(aes(xmin = Lower, xmax = Upper), height = 0,position=position_dodgev(height=-0.5)) +
  theme_bw() + theme(text=element_text(family="serif")) + 
  ylab("") + xlab("Effect of Protester Identity (Republican=Right)") + ggtitle("")
dev.off()

# Figure 1: stacked bar graph of raw results

bar_stats <- data.frame(Answer=character(length = 20),Ans=numeric(length=20),
                        Protesters = character(length=20), Respondents=character(length=20), N =numeric(length=20))

bar_stats$Ans <- rep(1:5,times=4)
bar_stats$Answer <- rep(c("Strongly Disagree","Somewhat Disagree","Neither","Somewhat Agree","Strongly Agree"),times=4)
bar_stats$Protesters <- rep(c("Democrats","Republicans","Democrats","Republicans"),each=5)
bar_stats$Respondents <- rep(c("Dem Respondents","GOP Respondents"),each=10)


for(i in 1:20){
  typegroup <- comp_dat[comp_dat$military==0&comp_dat$treat_1==ifelse(bar_stats$Protesters[i]=="Republicans",1,0)& 
                     comp_dat$pid_lean==ifelse(bar_stats$Respondents[i]=="GOP Respondents",1,0),]
  bar_stats$N[i] <- mean(typegroup$vig1_army==bar_stats$Ans[i],na.rm=T)*100
}

bar_stats$Answer <- factor(bar_stats$Answer, levels=c("Strongly Disagree","Somewhat Disagree","Neither","Somewhat Agree","Strongly Agree"))
bar_stats$Protesters <- factor(bar_stats$Protesters,levels=c("Democrats","Republicans"))
bar_stats$Respondents <- factor(bar_stats$Respondents,levels=c("Dem Respondents","GOP Respondents"))

pdf("Figures/Figure_1a.pdf", width = 5, height = 2.5)
ggplot(data=bar_stats, aes(x=Protesters, y=N,group=Answer, color=Answer,fill=Answer)) +
  ylab("") + xlab("Party of Protesters")+  labs(color="Deploy Army?",fill="Deploy Army?")+
  theme(text=element_text(family="serif")) +  facet_wrap( ~ Respondents, nrow=2 ) +
  scale_fill_brewer(palette = "RdYlBu") + scale_colour_brewer(palette = "RdYlBu") +
  geom_bar(stat="identity") + scale_y_continuous(limits = c(0, 100)) + coord_flip()
dev.off()

# Same graph for Vignette 2

for(i in 1:20){
  typegroup <- comp_dat[comp_dat$military==0& comp_dat$treat_2==ifelse(bar_stats$Protesters[i]=="Republicans",1,0)& 
                          comp_dat$pid_lean==ifelse(bar_stats$Respondents[i]=="GOP Respondents",1,0),]
  bar_stats$N[i] <- mean(typegroup$vig2_army==bar_stats$Ans[i],na.rm=T)*100
}

pdf("Figures/Figure_1b.pdf", width = 5, height = 2.5)
ggplot(data=bar_stats, aes(x=Protesters, y=N,group=Answer, color=Answer,fill=Answer)) +
  ylab("") + xlab("Party of Protesters")+  labs(color="Deploy Army?",fill="Deploy Army?")+
  theme(text=element_text(family="serif")) +  facet_wrap( ~ Respondents,nrow=2 ) +
  scale_fill_brewer(palette = "RdYlBu") + scale_colour_brewer(palette = "RdYlBu") +
  geom_bar(stat="identity") + scale_y_continuous(limits = c(0, 100)) + coord_flip()
dev.off()


# Figure A1: Raw Preference for Intervention by PID

results <- data.frame(Vignette=numeric(length=18),Vignette_Name=character(length=18),DV=character(length=18),DV_Name=character(length=18),
                      Resp=numeric(length=18),Respondents=character(length=18),Mean=numeric(length=18),SE=numeric(length=18))

results$Vignette <- rep(1:2, each=9)
results$Vignette_Name <- rep(c("Capitol Protest","Downtown Riot"),each=9)
results$Resp <- rep(1:3,times=6)
results$Respondents <- rep(c("Democrats","Independents","Republicans"),times=6)
results$DV_Name <- rep(rep(c("Deploy Army","Militarize Police","Armed Volunteers"),each=3),times=2)
results$DV <- rep(1:6,each=3)

all_vars <- cbind(comp_dat$vig1_army>=4,ifelse(comp_dat$vig1_police_eq%in%c(1,3),1,0),comp_dat$vig1_militia>=4,
                  comp_dat$vig2_army>=4,ifelse(comp_dat$vig2_police_eq%in%c(1,3),1,0),comp_dat$vig2_militia>=4)

se <- function(x) sqrt(var(x)/length(x))

for(i in 1:18){
  dv <- all_vars[comp_dat$military==0&comp_dat$pid_bin==results$Resp[i],results$DV[i]]
  dv <- dv[!is.na(dv)]*100
  results[i,7] <- mean(dv)
  results[i,8] <- se(dv)
}

results$Upper <- results$Mean+1.96*results$SE
results$Lower <- results$Mean-1.96*results$SE
results$Value <- as.character(round(results$Mean))

results$Respondents <- factor(results$Respondents,levels=c("Republicans","Independents","Democrats"))
results$DV_Name <- factor(results$DV_Name,levels=c("Armed Volunteers","Militarize Police","Deploy Army"))


pdf("Figures/Figure_A1.pdf", width = 7, height = 4)
ggplot(results, aes(x = Mean, y = DV_Name, color=Respondents,group=Respondents,label=Value)) +
  facet_grid(Vignette_Name ~.) +
  geom_point(size = 2,position=position_dodgev(height =-0.5)) + 
  scale_x_continuous(limits = c(0, 100)) +
  geom_errorbar(aes(xmin = Lower, xmax = Upper), width=0,position=position_dodgev(height =-0.5)) +
  theme_bw() + theme(text=element_text(family="serif")) + 
  labs(color="Respondents",color="Respondents")+ geom_text(family="serif",hjust=-2, vjust=.5,position=position_dodgev(height =-.5)) +
  ylab(" ") + xlab("Percent Who Support") + ggtitle("")
dev.off()


# Figure 4: graph of raw preferences by party (binary)

results <- data.frame(Vignette=numeric(length=12),Vignette_Name=character(length=12),DV=character(length=12),DV_Name=character(length=12),
                      Resp=numeric(length=12),Respondents=character(length=12),Mean=numeric(length=12),SE=numeric(length=12))

results$Vignette <- rep(1:2, each=6)
results$Vignette_Name <- rep(c("Capitol Protest","Downtown Riot"),each=6)
results$Resp <- rep(1:2,times=6)
results$Respondents <- rep(c("Democrats","Republicans"),times=6)
results$DV_Name <- rep(rep(c("Deploy Army","Militarize Police","Armed Volunteers"),each=2),times=2)
results$DV <- rep(1:6,each=2)

all_vars <- cbind(comp_dat$vig1_army>=4,ifelse(comp_dat$vig1_police_eq%in%c(1,3),1,0),comp_dat$vig1_militia>=4,
                  comp_dat$vig2_army>=4,ifelse(comp_dat$vig2_police_eq%in%c(1,3),1,0),comp_dat$vig2_militia>=4)

for(i in 1:12){
  dv <- all_vars[comp_dat$military==0&comp_dat$pid_lean==ifelse(results$Respondents[i]=="Democrats",0,1),results$DV[i]]
  dv <- dv[!is.na(dv)]*100
  results[i,7] <- mean(dv)
  results[i,8] <- se(dv)
}

results$Upper <- results$Mean+1.96*results$SE
results$Lower <- results$Mean-1.96*results$SE
results$Value <- as.character(round(results$Mean))

results$Respondents <- factor(results$Respondents,levels=c("Republicans","Democrats"))
results$DV_Name <- factor(results$DV_Name,levels=c("Armed Volunteers","Militarize Police","Deploy Army"))


pdf("Figures/Figure_4.pdf", width = 7, height = 3)
ggplot(results, aes(x = Mean, y = DV_Name, color=Respondents,group=Respondents,label=Value)) +
  facet_grid(Vignette_Name ~.) +
  geom_point(size = 2,position=position_dodgev(height =-.5)) + 
  scale_x_continuous(limits = c(0, 100)) +
  geom_errorbar(aes(xmin = Lower, xmax = Upper), width=0,position=position_dodgev(height =-.5)) +
  theme_bw() + theme(text=element_text(family="serif")) + 
  labs(color="Respondents",color="Respondents")+ geom_text(family="serif",hjust=-2, vjust=.5,position=position_dodgev(height =-.5)) +
  ylab("") + xlab("Percent Who Support") + ggtitle("")
dev.off()


# Figure 2: Trust in Institutions Figure

bar_stats <- data.frame(Answer=character(length = 24),Ans=numeric(length=24),
                        Institution = character(length=24), Respondents=character(length=24), N =numeric(length=24))

bar_stats$Ans <- rep(1:4,times=6)
bar_stats$Answer <- rep(c("Just about never","Some of the time","Most of the time","Just about always"),times=6)
bar_stats$Institution <- rep(c("Fed Government","US Military","Police"),each=8)
bar_stats$Inst <- rep(1:3,each=8)
bar_stats$Respondents <- rep(rep(c("Dem","Rep"),each=4),times=3)

trust_vars <- cbind(comp_dat$trust_prez,comp_dat$trust_military,comp_dat$trust_police)

for(i in 1:24){
  typegroup <- trust_vars[comp_dat$military==0 & !is.na(comp_dat$pid_lean) &
                            comp_dat$pid_lean==ifelse(bar_stats$Respondents[i]=="Rep",1,0),
                          bar_stats$Inst[i]]
  bar_stats$N[i] <- mean(typegroup==bar_stats$Ans[i],na.rm=T)*100
}

bar_stats$Answer <- factor(bar_stats$Answer, levels=c("Just about never","Some of the time","Most of the time","Just about always"))
bar_stats$Institution <- factor(bar_stats$Institution,levels=c("Fed Government","US Military","Police"))
bar_stats$Respondents <- factor(bar_stats$Respondents,levels=c("Dem","Rep"))

pdf("Figures/Figure_2.pdf", width = 6, height = 3)
ggplot(data=bar_stats, aes(x=Respondents, y=N,group=Answer, color=Answer,fill=Answer)) +
  ylab("") + xlab("Respondent Party ID (+ Leaners)")+  labs(color="How often can you trust...?",fill="How often can you trust...?")+
  theme(text=element_text(family="serif")) +  facet_wrap( ~ Institution,ncol=1) +
  scale_fill_brewer(palette = "RdYlBu") + scale_colour_brewer(palette = "RdYlBu") +
  geom_bar(stat="identity") + scale_y_continuous(limits = c(0, 100)) + coord_flip()
dev.off()

# Table A1: Trust in Institutions

trust_vars <- cbind(comp_dat$trust_prez,comp_dat$trust_congress,comp_dat$trust_military,comp_dat$trust_police)
groups <- ifelse(comp_dat$military==1,4,comp_dat$pid_bin)
trust_mat <- matrix(NA,nrow=4,ncol=4)
colnames(trust_mat) <- c("Democrats","Independents","Republicans","Military")
rownames(trust_mat) <- c("Federal Govt","State Govt","Military","Local Police")
for(i in 1:4){
  for(j in 1:4){
    trust_mat[i,j] <- mean(trust_vars[groups==j,i]>=3,na.rm=T)
  }
}
round(trust_mat*100,0)

# Figure A3: correlation between trust and intervention support

comp_dat$copartisan_1 <- ifelse(comp_dat$treat_1==comp_dat$pid_lean,1,0)
comp_dat$copartisan_2 <- ifelse(comp_dat$treat_2==comp_dat$pid_lean,1,0)
comp_dat$trust_simple <- (comp_dat$trust_military==4)*3+(comp_dat$trust_military==3)*2+(comp_dat$trust_military%in%1:2)*1

results <- data.frame(Vignette=numeric(length=12),Vignette_Name=character(length=12),Trust=character(length=12),Resp=numeric(length=12),
                      Treat=numeric(length=12),Treatment=character(length=12),Mean=numeric(length=12),SE=numeric(length=12))

results$Vignette <- rep(1:2, each=6)
results$Vignette_Name <- rep(c("Capitol Protest","Downtown Riot"),each=6)
results$Resp <- rep(rep(1:3,each=2),times=2)
results$Trust <- rep(rep(c("Low","Medium","High"),each=2),times=2)
results$Treat <- rep(0:1,times=6)
results$Treatment <- rep(c("Opposite Party","Same Party"),times=6)

all_vars <- cbind(comp_dat$vig1_army>=4,comp_dat$vig2_army>=4)
all_copartisan <- cbind(comp_dat$same_1,comp_dat$same_2)

se <- function(x) sqrt(var(x)/length(x))

for(i in 1:12){
  dv <- all_vars[comp_dat$military==0&comp_dat$trust_simple==results$Resp[i]&all_copartisan[,results$Vignette[i]]==results$Treat[i],
                 results$Vignette[i]]
  dv <- dv[!is.na(dv)]*100
  results[i,7] <- mean(dv)
  results[i,8] <- se(dv)
}

results$Upper <- pmin(results$Mean+1.96*results$SE,100)
results$Lower <- pmax(results$Mean-1.96*results$SE,0)

results$Vignette_Name <- factor(results$Vignette_Name,levels = c("Capitol Protest","Downtown Riot"))
results$Trust <- factor(results$Trust,levels=c("Low","Medium","High"))
results$Treatment <- factor(results$Treatment,levels=c("Same Party","Opposite Party"))
results$Value <- as.character(round(results$Mean))

pdf("Figures/Figure_A3.pdf", width = 7, height = 3)
ggplot(results, aes(x = Mean, y = Trust, color=Treatment,group=Treatment,label=Value)) +
  facet_grid(Vignette_Name ~.) +
  geom_point(size = 2,position=position_dodgev(height =-0.5)) + 
  scale_x_continuous(limits = c(0, 100)) +
  geom_errorbar(aes(xmin = Lower, xmax = Upper), width=0,position=position_dodgev(height =-0.5)) +
  theme_bw() + theme(text=element_text(family="serif")) + 
  labs(color="Protesters",color="Protesters")+ geom_text(family="serif",hjust=-3, vjust=.5,position=position_dodgev(height =-.5)) +
  ylab("Trust in Military") + xlab("Percent Who Support") + ggtitle("")
dev.off()


# Figure 6: Military vs Public Support for Interventions

results <- data.frame(Vignette=numeric(length=12),Vignette_Name=character(length=12),DV=character(length=12),DV_Name=character(length=12),
                      Resp=numeric(length=12),Respondents=character(length=12),Mean=numeric(length=12),SE=numeric(length=12))

results$Vignette <- rep(1:2, each=6)
results$Vignette_Name <- rep(c("Capitol Protest","Downtown Riot"),each=6)
results$Resp <- rep(1:2,times=6)
results$Respondents <- rep(c("Military","Public"),times=6)
results$DV_Name <- rep(rep(c("Deploy Army","Militarize Police","Armed Volunteers"),each=2),times=2)
results$DV <- rep(1:6,each=2)

all_vars <- cbind(comp_dat$vig1_army>=4,ifelse(comp_dat$vig1_police_eq%in%c(1,3),1,0),comp_dat$vig1_militia>=4,
                  comp_dat$vig2_army>=4,ifelse(comp_dat$vig2_police_eq%in%c(1,3),1,0),comp_dat$vig2_militia>=4)

se <- function(x) sqrt(var(x)/length(x))

for(i in 1:12){
  dv <- all_vars[comp_dat$military==(2-results$Resp[i]),results$DV[i]]
  dv <- dv[!is.na(dv)]*100
  results[i,7] <- mean(dv)
  results[i,8] <- se(dv)
}

results$Upper <- pmin(results$Mean+1.96*results$SE,100)
results$Lower <- pmax(results$Mean-1.96*results$SE,0)

results$Respondents <- factor(results$Respondents,levels=c("Military","Public"))
results$DV_Name <- factor(results$DV_Name,levels=c("Armed Volunteers","Militarize Police","Deploy Army"))
results$Value <- as.character(round(results$Mean))

pdf("Figures/Figure_6.pdf", width = 7, height = 3)
ggplot(results, aes(x = Mean, y = DV_Name, color=Respondents,group=Respondents,label=Value)) +
  facet_grid(Vignette_Name ~.) +
  geom_point(size = 2,position=position_dodgev(height =-0.5)) + 
  scale_x_continuous(limits = c(0, 100)) +
  geom_errorbar(aes(xmin = Lower, xmax = Upper), width=0,position=position_dodgev(height =-0.5)) +
  theme_bw() + theme(text=element_text(family="serif")) + 
  labs(color="Respondents",color="Respondents")+ geom_text(family="serif",hjust=-3, vjust=.5,position=position_dodgev(height =-.5)) +
  ylab("") + xlab("Percent Who Support") + ggtitle("")
dev.off()


# Figure 5: Raw Support for Army Deployment, Mil vs Civilian

all_vars <- cbind(comp_dat$vig1_army,comp_dat$vig2_army)

bar_stats <- data.frame(Answer=character(length = 20),Ans=numeric(length=20), Q=numeric(length=20),
                        Question = character(length=20), Respondents=character(length=20), N =numeric(length=20))

bar_stats$Ans <- rep(1:5,times=4)
bar_stats$Answer <- rep(c("Strongly Disagree","Somewhat Disagree","Neither","Somewhat Agree","Strongly Agree"),times=4)
bar_stats$Question <- rep(c("Capitol Protest","Downtown Riots"),each=10)
bar_stats$Q <- rep(1:2,each=10)
bar_stats$Respondents <- rep(rep(c("Public","Military"),each=5),times=2)

for(i in 1:20){
  typegroup <- all_vars[comp_dat$military==ifelse(bar_stats$Respondents[i]=="Military",1,0),
                        bar_stats$Q[i]]
  bar_stats$N[i] <- mean(typegroup==bar_stats$Ans[i],na.rm=T)*100
}


bar_stats$Answer <- factor(bar_stats$Answer, levels=c("Strongly Disagree","Somewhat Disagree","Neither","Somewhat Agree","Strongly Agree"))
bar_stats$Question <- factor(bar_stats$Question,levels=c("Capitol Protest","Downtown Riots"))
bar_stats$Respondents <- factor(bar_stats$Respondents,levels=c("Public","Military"))

pdf("Figures/Figure_5.pdf", width = 5, height = 3)
ggplot(data=bar_stats, aes(x=Respondents, y=N,group=Answer, color=Answer,fill=Answer)) +
  ylab("") + xlab("Respondents")+  labs(color="Deploy Army?",fill="Deploy Army?")+
  theme(text=element_text(family="serif")) +  facet_wrap( ~ Question,nrow=2 ) +
  scale_fill_brewer(palette = "RdYlBu") + scale_colour_brewer(palette = "RdYlBu") +
  geom_bar(stat="identity") + scale_y_continuous(limits = c(0, 100))+coord_flip()
dev.off()


# Figure 3: Summary of Partisanship Effect on Civilians

#summary stats
outparty <- c(civ$pid_lean!=civ$treat_1,civ$pid_lean!=civ$treat_2)
army <- c(civ$vig1_army,civ$vig2_army)
police <- c(civ$vig1_police_eq%in%c(1,3),civ$vig2_police_eq%in%c(1,3))
militia <- c(civ$vig1_militia,civ$vig2_militia)

by(army>=4,INDICES=outparty,FUN=mean)
by(police,INDICES=outparty,FUN=mean)
by(militia>=4,INDICES=outparty,FUN=mean)


# Figure
results <- data.frame(Vignette=numeric(length=12),Vignette_Name=character(length=12),DV=character(length=12),DV_Name=character(length=12),
                      Resp=numeric(length=12),Respondents=character(length=12),Mean=numeric(length=12),SE=numeric(length=12))

results$Vignette <- rep(1:2, each=6)
results$Vignette_Name <- rep(c("Capitol Protest","Downtown Riot"),each=6)
results$Resp <- rep(1:2,times=6)
results$Respondents <- rep(c("Same Party","Opposite Party"),times=6)
results$DV_Name <- rep(rep(c("Deploy Army","Militarize Police","Armed Volunteers"),each=2),times=2)
results$DV <- rep(1:6,each=2)

all_vars <- cbind(comp_dat$vig1_army>=4,ifelse(comp_dat$vig1_police_eq%in%c(1,3),1,0),comp_dat$vig1_militia>=4,
                  comp_dat$vig2_army>=4,ifelse(comp_dat$vig2_police_eq%in%c(1,3),1,0),comp_dat$vig2_militia>=4)

for(i in 1:12){
  treat <- all_treats[comp_dat$military==0,results$Vignette[i]]
  party <- comp_dat$pid_lean[comp_dat$military==0]
  dv <- all_vars[comp_dat$military==0,results$DV[i]]
  sameparty <- ifelse(!is.na(party),treat==party,NA)
  dv <- dv[sameparty==(results$Respondents[i]=="Same Party")]
  dv <- dv[!is.na(dv)]*100
  results[i,7] <- mean(dv)
  results[i,8] <- se(dv)
}

results$Upper <- results$Mean+1.96*results$SE
results$Lower <- results$Mean-1.96*results$SE
results$Value <- as.character(round(results$Mean))

results$Respondents <- factor(results$Respondents,levels=c("Same Party","Opposite Party"))
results$DV_Name <- factor(results$DV_Name,levels=c("Armed Volunteers","Militarize Police","Deploy Army"))


pdf("Figures/Figure_3.pdf", width = 7, height = 3)
ggplot(results, aes(x = Mean, y = DV_Name, color=Respondents,group=Respondents,label=Value)) +
  facet_grid(Vignette_Name ~.) +
  geom_point(size = 2,position=position_dodgev(height =-.5)) + 
  scale_x_continuous(limits = c(0, 100)) +
  geom_errorbar(aes(xmin = Lower, xmax = Upper), width=0,position=position_dodgev(height =-.5)) +
  theme_bw() + theme(text=element_text(family="serif")) + 
  labs(color="Protester Party",color="Protester Party")+ geom_text(family="serif",hjust=-2, vjust=.5,position=position_dodgev(height =-.5)) +
  ylab("") + xlab("Percent Who Support") + ggtitle("")
dev.off()



# Figure A6: Partisan Treatment Effects on Military Personnel

results <- data.frame(Vignette=numeric(length=12),Vignette_Name=character(length=12),DV=character(length=12),DV_Name=character(length=12),
                      Resp=numeric(length=12),Respondents=character(length=12),Mean=numeric(length=12),SE=numeric(length=12))

results$Vignette <- rep(1:2, each=6)
results$Vignette_Name <- rep(c("Capitol Protest","Downtown Riot"),each=6)
results$Resp <- rep(1:2,times=6)
results$Respondents <- rep(c("Same Party","Opposite Party"),times=6)
results$DV_Name <- rep(rep(c("Deploy Army","Militarize Police","Armed Volunteers"),each=2),times=2)
results$DV <- rep(1:6,each=2)

all_vars <- cbind(comp_dat$vig1_army>=4,ifelse(comp_dat$vig1_police_eq%in%c(1,3),1,0),comp_dat$vig1_militia>=4,
                  comp_dat$vig2_army>=4,ifelse(comp_dat$vig2_police_eq%in%c(1,3),1,0),comp_dat$vig2_militia>=4)

for(i in 1:12){
  treat <- all_treats[comp_dat$military==1,results$Vignette[i]]
  party <- comp_dat$pid_lean[comp_dat$military==1]
  dv <- all_vars[comp_dat$military==1,results$DV[i]]
  sameparty <- ifelse(!is.na(party),treat==party,NA)
  dv <- dv[sameparty==(results$Respondents[i]=="Same Party")]
  dv <- dv[!is.na(dv)]*100
  results[i,7] <- mean(dv)
  results[i,8] <- se(dv)
}

results$Upper <- pmin(results$Mean+1.96*results$SE,100)
results$Lower <- pmax(results$Mean-1.96*results$SE,0)
results$Value <- as.character(round(results$Mean))

results$Respondents <- factor(results$Respondents,levels=c("Same Party","Opposite Party"))
results$DV_Name <- factor(results$DV_Name,levels=c("Armed Volunteers","Militarize Police","Deploy Army"))


pdf("Figures/Figure_A6.pdf", width = 7, height = 3)
ggplot(results, aes(x = Mean, y = DV_Name, color=Respondents,group=Respondents,label=Value)) +
  facet_grid(Vignette_Name ~.) +
  geom_point(size = 2,position=position_dodgev(height =-.5)) + 
  scale_x_continuous(limits = c(0, 100)) +
  geom_errorbar(aes(xmin = Lower, xmax = Upper), width=0,position=position_dodgev(height =-.5)) +
  theme_bw() + theme(text=element_text(family="serif")) + 
  labs(color="Protester Party",color="Protester Party")+ geom_text(family="serif",hjust=-2, vjust=.5,position=position_dodgev(height =-.5)) +
  ylab("") + xlab("Percent Who Support") + ggtitle("")
dev.off()


# Figure A5: Partisanship effects on strong partisans

results <- data.frame(Vignette=numeric(length=12),Vignette_Name=character(length=12),DV=character(length=12),DV_Name=character(length=12),
                      Resp=numeric(length=12),Respondents=character(length=12),Mean=numeric(length=12),SE=numeric(length=12))

results$Vignette <- rep(1:2, each=6)
results$Vignette_Name <- rep(c("Capitol Protest","Downtown Riot"),each=6)
results$Resp <- rep(1:2,times=6)
results$Respondents <- rep(c("Same Party","Opposite Party"),times=6)
results$DV_Name <- rep(rep(c("Deploy Army","Militarize Police","Armed Volunteers"),each=2),times=2)
results$DV <- rep(1:6,each=2)

all_vars <- cbind(comp_dat$vig1_army>=4,ifelse(comp_dat$vig1_police_eq%in%c(1,3),1,0),comp_dat$vig1_militia>=4,
                  comp_dat$vig2_army>=4,ifelse(comp_dat$vig2_police_eq%in%c(1,3),1,0),comp_dat$vig2_militia>=4)

for(i in 1:12){
  treat <- all_treats[comp_dat$military==0&comp_dat$pid_scale%in%c(1,7),results$Vignette[i]]
  party <- comp_dat$pid_lean[comp_dat$military==0&comp_dat$pid_scale%in%c(1,7)]
  dv <- all_vars[comp_dat$military==0&comp_dat$pid_scale%in%c(1,7),results$DV[i]]
  sameparty <- ifelse(!is.na(party),treat==party,NA)
  dv <- dv[sameparty==(results$Respondents[i]=="Same Party")]
  dv <- dv[!is.na(dv)]*100
  results[i,7] <- mean(dv)
  results[i,8] <- se(dv)
}

results$Upper <- results$Mean+1.96*results$SE
results$Lower <- results$Mean-1.96*results$SE
results$Value <- as.character(round(results$Mean))

results$Respondents <- factor(results$Respondents,levels=c("Same Party","Opposite Party"))
results$DV_Name <- factor(results$DV_Name,levels=c("Armed Volunteers","Militarize Police","Deploy Army"))


pdf("Figures/Figure_A5.pdf", width = 7, height = 3)
ggplot(results, aes(x = Mean, y = DV_Name, color=Respondents,group=Respondents,label=Value)) +
  facet_grid(Vignette_Name ~.) +
  geom_point(size = 2,position=position_dodgev(height =-.5)) + 
  scale_x_continuous(limits = c(0, 100)) +
  geom_errorbar(aes(xmin = Lower, xmax = Upper), width=0,position=position_dodgev(height =-.5)) +
  theme_bw() + theme(text=element_text(family="serif")) + 
  labs(color="Protester Party",color="Protester Party")+ geom_text(family="serif",hjust=-2, vjust=.5,position=position_dodgev(height =-.5)) +
  ylab("") + xlab("Percent Who Support") + ggtitle("")
dev.off()


# Figure 7: Mechanisms of Effects on Civilians

groups <- cbind(comp_dat$military, 1-comp_dat$military,
                ifelse(!is.na(comp_dat$mil_connections)&comp_dat$mil_connections==1,1,0),
                ifelse(!is.na(comp_dat$vet_status)& comp_dat$vet_status==2,1,0),
                ifelse(!is.na(comp_dat$gender)&comp_dat$age>=25&comp_dat$age<=55&
                         comp_dat$gender==1&comp_dat$education>=5,1,0))

results <- data.frame(Resp=numeric(length=ncol(groups)),Respondents=character(length=ncol(groups)),
                      N=numeric(length=ncol(groups)),
                      Mean=numeric(length=ncol(groups)),SE=numeric(length=ncol(groups)))

results$Resp <- 1:5
results$Respondents <- c("Military Officers","Public","Friends & Family","Veterans","Men w/ BA")
results$Group <- c("Military","Public","Subgroup","Subgroup","Subgroup")
results$Group <- factor(results$Group,levels=c("Military","Subgroup","Public"))

for(i in 1:5){
  grp <- groups[,i]
  dv <- c(ifelse(comp_dat$vig1_army[grp==1]>=4,1,0),ifelse(comp_dat$vig2_army[grp==1]>=4,1,0))*100
  dv <- dv[!is.na(dv)]
  results$N[i] <- length(dv)/2
  results$Mean[i] <- mean(dv,na.rm = T)
  results$SE[i] <- se(dv)
}

results$Upper <- pmin(results$Mean+1.96*results$SE,100)
results$Lower <- pmax(results$Mean-1.96*results$SE,0)

results$Respondents <- factor(results$Respondents,
                              levels=c("Men w/ BA","Friends & Family","Veterans","Public","Military Officers"))

baseline1 <- results$Mean[1]
baseline2 <- results$Mean[2]

pdf("Figures/Figure_7.pdf", width = 4.5, height = 2)
ggplot(results, aes(x = Mean, y = Respondents,color=Group,fill=Group)) +
  geom_point(size=2) + 
  scale_x_continuous(limits = c(0, 100)) +
  geom_errorbar(aes(xmin = Lower, xmax = Upper),width=0) + 
  geom_vline(xintercept = baseline1,color="red",linetype="dashed") + geom_vline(xintercept = baseline2,color="blue",linetype="dashed") +
  theme_bw() + theme(text=element_text(family="serif"),legend.position="none") + 
  ylab("Respondents") + xlab("Percent Who Support Deploying the Army") + ggtitle("")
dev.off()


# Figure A7: graph of support by tenure in military

comp_dat$older <- comp_dat$tenure>2

results <- data.frame(Vignette=numeric(length=12),Vignette_Name=character(length=12),DV=character(length=12),DV_Name=character(length=12),
                      Resp=numeric(length=12),Respondents=character(length=12),Mean=numeric(length=12),SE=numeric(length=12))

results$Vignette <- rep(1:2, each=6)
results$Vignette_Name <- rep(c("Capitol Protest","Downtown Riot"),each=6)
results$Resp <- rep(1:2,times=6)
results$Respondents <- rep(c("Shorter-Tenured","Longer-Tenured"),times=6)
results$DV_Name <- rep(rep(c("Deploy Army","Militarize Police","Armed Volunteers"),each=2),times=2)
results$DV <- rep(1:6,each=2)

all_vars <- cbind(comp_dat$vig1_army>=4,ifelse(comp_dat$vig1_police_eq%in%c(1,3),1,0),comp_dat$vig1_militia>=4,
                  comp_dat$vig2_army>=4,ifelse(comp_dat$vig2_police_eq%in%c(1,3),1,0),comp_dat$vig2_militia>=4)

for(i in 1:12){
  dv <- all_vars[comp_dat$military==1&comp_dat$older==ifelse(results$Respondents[i]=="Longer-Tenured",1,0),results$DV[i]]
  dv <- dv[!is.na(dv)]*100
  results[i,7] <- mean(dv)
  results[i,8] <- se(dv)
}

results$Upper <- results$Mean+1.96*results$SE
results$Lower <- pmax(results$Mean-1.96*results$SE,rep(0,12))
results$Value <- as.character(round(results$Mean))

results$Respondents <- factor(results$Respondents,levels=c("Longer-Tenured","Shorter-Tenured"))
results$DV_Name <- factor(results$DV_Name,levels=c("Armed Volunteers","Militarize Police","Deploy Army"))


pdf("Figures/Figure_A7.pdf", width = 7, height = 3)
ggplot(results, aes(x = Mean, y = DV_Name, color=Respondents,group=Respondents,label=Value)) +
  facet_grid(Vignette_Name ~.) +
  geom_point(size = 2,position=position_dodgev(height =-.5)) + 
  scale_x_continuous(limits = c(0, 100)) +
  geom_errorbar(aes(xmin = Lower, xmax = Upper), width=0,position=position_dodgev(height =-.5)) +
  theme_bw() + theme(text=element_text(family="serif")) + 
  labs(color="Respondents",color="Respondents")+ geom_text(family="serif",hjust=-2, vjust=.5,position=position_dodgev(height =-.5)) +
  ylab("") + xlab("Percent Who Support") + ggtitle("")
dev.off()


## Figure A4: Attention Check

comp_dat$attention_1 <- ifelse(comp_dat$treat_1==0,comp_dat$vig1_gop_atten_check,comp_dat$vig1_dem_atten_check)
comp_dat$attention_2 <- ifelse(comp_dat$treat_2==1,comp_dat$vig2_gop_atten_check,comp_dat$vig2_dem_atten_check)
table(comp_dat$treat_1,comp_dat$attention_1)
table(comp_dat$treat_2,comp_dat$attention_2)

bar_stats <- data.frame(Answer=character(length = 8),Ans=numeric(length=8),
                        Treatment = character(length=8), Respondents=character(length=8), N =numeric(length=8))

bar_stats$Ans <- rep(1:2,times=4)
bar_stats$Answer <- rep(c("Republican","Democrat"),times=4)
bar_stats$Treatment <- rep(c("Democrat","Democrat","Republican","Republican"),times=2)
bar_stats$Respondents <- rep(c("Civilian Respondents","Military Respondents"),each=4)

for(i in 1:8){
  typegroup <- comp_dat$attention_1[comp_dat$military==ifelse(bar_stats$Respondents[i]=="Military Respondents",1,0)&comp_dat$treat_1==ifelse(bar_stats$Treatment[i]=="Democrat",1,0)]
  bar_stats$N[i] <- mean(typegroup==bar_stats$Ans[i],na.rm=T)*100
}

bar_stats$Answer <- factor(bar_stats$Answer, levels=c("Republican","Democrat"))
bar_stats$Treatment <- factor(bar_stats$Treatment,levels=c("Republican","Democrat"))
bar_stats$Respondents <- factor(bar_stats$Respondents,levels=c("Civilian Respondents","Military Respondents"))
bar_stats$N <- round(bar_stats$N)
bar_stats

pdf("Figures/Figure_A4.pdf", width = 5, height = 3)
ggplot(data=bar_stats, aes(x=Treatment, y=N,group=Answer, color=Answer,fill=Answer)) +
  ylab("Given Answer (Election Winner)") + xlab("Correct Answer")+  labs(color="",fill="")+
  theme(text=element_text(family="serif")) +  facet_wrap( ~ Respondents, nrow=2 ) +
  scale_fill_brewer(palette = "Set1") + scale_colour_brewer(palette = "Set1") +
  geom_bar(stat="identity") + scale_y_continuous(limits = c(0, 100)) + coord_flip()
dev.off()

